The time-reversal test for stochastic quantum dynamics 
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The calculation of quantum dynamics is currently a central issue in theoretical physics, with diverse appli- 
cations ranging from ultra-cold atomic Bose-Einstein condensates (BEC) to condensed matter, biology, and 
even astrophysics. Here we demonstrate a conceptually simple method of determining the regime of validity 
of stochastic simulations of unitary quantum dynamics by employing a time-reversal test. We apply this test to 
a simulation of the evolution of a quantum anharmonic oscillator with up to 6.022 x 10 23 (Avogadro's num- 
ber) of particles. This system is realisable as a Bose-Einstein condensate in an optical lattice, for which the 
time-reversal procedure could be implemented experimentally. 



The difficulty of real-time quantum dynamical calculations 
is caused by complexity. The computational resources re- 
quired to directly represent the Hilbert space of a large quan- 
tum system are enormous. This problem led to Feynman's 
proposal 1 1 ] to develop quantum computer hardware for quan- 
tum dynamics. In the absence of such devices, digital comput- 
ers must be employed for these calculations at present. Re- 
gardless of the hardware or software being utilised, there is a 
profound question of how one can check that results are calcu- 
lated accurately. This is an especially difficult issue with the 
time-evolution of quantum many-body systems, which is one 
of the central challenges in theoretical physics. There are few 
exact solutions, yet results must be calculated systematically 
to within a known error in order to allow theoretical predic- 
tions to be tested experimentally. 

Stimulated by the success of quantum Monte-Carlo meth- 
ods in imaginary time B BBj the method used here for 
real-time quantum dynamics relies on sampling a probabilis- 
tic phase-space representation. Related approaches include 
Wigner's classical phase-space representation |6], which was 
used to develop semi-classical approximations similar to those 
for quantum chaos calculations 17], as well as other classi- 
cal phase-space B HUH Eh representations. More recent 
phase-space methods for quantum simulations use a nonclas- 
sical phase-space together with a weight parameter analogous 
to those used in path-integrals II 211 . These methods allow 
quantum dynamical simulations from first principles without 
semiclassical approximations. However, the sampling error 
can become a limiting factor. 

Fortunately, an important property of time-independent 
Hamiltonians is that evolution backward in time is equivalent 
to evolution forward in time under a Hamiltonian of the op- 
posite sign. This suggests a simple yet powerful test that any 
unitary quantum dynamical simulation must pass. Beginning 
with a well-defined initial state, a simulation is evolved for a 
time period for which we are interested in the quantum dy- 
namics. The Hamiltonian is then negated and the simulation 
evolved again for the same period. For reliable simulation all 
relevant initial observables should be recovered. 

Phase-space methods utilising quasi-probability distribu- 
tions lead one to sample an equivalent set of stochastic differ- 
ential equations (SDEs) with random noise terms, and these 



techniques scale linearly with the number of modes 

11611 . While such methods have been successful for many 
problems the sampling errors sometimes grow in 

time and eventually can become unmanageable. Similar is- 
sues are encountered in simulating classical chaos, where sen- 
sitive dependence on initial conditions leading to an expo- 
nential growth of errors fl9ll can be tested via time-reversal. 
However, the use of intrinsically random equations for time- 
reversible quantum evolution appears paradoxical. How can 
one have time-reversibility in a method which appears to in- 
troduce increasing entropy at each step? It is this question we 
focus on here, by showing that this type of time-evolution is 
in fact completely reversible due to the storage of information 
in quantum correlations. 

All currently known phase-space methods can be repre- 
sented in a unified manner by an expansion of the density op- 
erator as 



daG(a)A(a), 



(1) 



where G(a) is a positive distribution function over the phase- 
space a, and A(a) is an over-complete basis for the Hilbert 
space 1 20]. A variety of techniques can be realised by 
changing the basis set, the dynamical equations (which are 
equivalent under a "stochastic gauge" symmetry |21]), and 
the numerical integration algorithm. We illustrate the time- 
reversal test for the particular case of a stochastic gauge 
simulation Elll . For this method the phase-space is 
<3 = 0r f2), which is a 2M + 1 complex dimensional 
vector containing phase-space variables a k and (3k (where 
a = {a±, . . . , ctk, ■ ■ ■ , olm}, etc.) for each of M bosonic 
modes, together with an additional variable Q, termed the 
weight. The operator basis is 

W){Pl\ 



A(a) = Q( 



(2) 



where the coherent state ja*.) is an eigenstate of the boson 
annihilation operator a/, for the kth mode, with a mean boson 
number fife = (a\.a,k) = \ctk\ 2 - 

For two-body interactions the master equation for time evo- 
lution of the density operator can be shown to be equivalent 
to a Fokker-Planck equation for the evolution of the quasi- 
probability gauge distribution G(a) with basis 0. This in 
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turn is equivalent to a set of SDEs. The moments of the gauge 
distribution function are then equivalent to dynamical quan- 
tum averages of products of bosonic creation and annihila- 
tion operators. For a single mode this equivalence can be ex- 
pressed (from now on we omit the mode label k) as 



<(« f ) 



;qm 



(Qf3 m a n 



stoch 



(n + n* 



1 stoch 



(3) 



where ()qm indicates a quantum mechanical average, and 
Ostoch is a stochastic average. 

In this paper we consider the anharmonic oscillator Hamil- 
tonian 



- Hk ^+2-2 

H = — a' a , 
2 



(4) 



which describes both the Kerr effect in nonlinear optics, and 
a single mode Bose-Einstein condensate (BEC). It is perhaps 
the simplest model of a many -body quantum system, and, as it 
is analytically solvable, it provides an excellent testing ground 
for simulation methods. 

An important quantum feature of this Hamiltonian is that 
given an initial coherent state p(Q) = \a) (a\, the dynam- 
ics display a series of collapses and revivals. From the an- 
alytic solution it is known that the quantum averages of the 



quadrature variables X 



at)/2 andF 



aJ)/2i 



undergo oscillations that eventually damp out to zero. How- 
ever, after a certain time the oscillations revive, and the ini- 
tial state is recovered. Defining the dimensionless time vari- 
able t = yfTiKijl-K, where n = \a\ 2 is the mean particle 
number, the relevant time scales are the oscillation period 
t sc ~ 0(l/\/E), the collapse time t co ii ~ 0(1), and the 
revival time r rev = v^- F° r large n the revival time is many 
times the collapse time, which in turn is much longer than the 
natural oscillation period. This revival is a uniquely quantum 
feature that does not occur in classical dynamics. 

While single mode Hamiltonians are often not a good de- 
scription of real systems, the anharmonic oscillator can be a 
good approximation for a Bose-Einstein condensate (BEC) in 
an optical lattice in the Mott insulator regime ll22il . A sud- 
den increase in the lattice depth from the superfluid regime 
can create coherent superpositions of atoms at each site which 
can be approximated by a coherent state. Indeed, such col- 
lapses and revivals have been observed with a BEC in a deep 
lattice lEl . 

With a suitable choice of stochastic gauge representa- 
tion 1 12| it was found to be possible to simulate past the col- 
lapse time with small statistical error using a modest number 
(~ 10 4 ) of stochastic trajectories 12111 . Here we check this 
calculation with the time-reversal test to demonstrate that the 
full quantum nature of the dynamics is preserved, even when 
the mean quadrature amplitudes are near zero. 

One finds |21] that the Ito SDEs corresponding to the an- 
harmonic oscillator Hamiltonian (0} are 



a = ia [-Kn x + £ 9! i] , 
/3 = i(3 [kji x + £* >2 ] , 



(5) 
(6) 




FIG. 1: (Color online) Time-reversal test for a stochastic gauge simu- 
lation of the X -quadrature of the anharmonic oscillator, (a) An initial 
coherent state with mean boson number n — 100 and 10° stochastic 
trajectories, (b) An initial coherent state with n = 6.022 x 10 23 and 
10 7 stochastic trajectories, with the calculation performed in the ro- 
tating frame. The solid lines are the simulation result for (A"), with 
the statistical error bars (representing one standard deviation) shown 
by the grey shading. For both cases the time-reversal was imple- 
mented at tr = 0.5 by negating the sign of the Hamiltonian and the 
initial state is recovered to within statistical error at r = 2tr. The 
forward time evolution to 2tr demonstrating the collapse to (X) = 
is also shown in both (a) and (b). 



in [£i - i£ 2 ] 



(7) 



where n = n x 



— a/3 is a complex variable correspond- 



ing to the particle number. Here £ g j are defined as transfor- 
mations of the fundamental noise terms £j through the intro- 
duction of an arbitrary stochastic diffusion gauge g, chosen 
for efficiency: 



(8) 



£ g ,j(t) = Vik [£j coshg + sxnhg] . 



For numerical integration with time steps dt, the £j can be 
implemented by independent real Gaussian noises of variance 
1 /dt and mean zero at each time step. The drift gauge used to 
obtain the deterministic parts of the equations from the Hamil- 
tonian has been described previously 12111 . It is convenient to 
transform these equations to logarithmic variables, and to use 
the Stratonovich calculus for integration |14]. To obtain the 
time-reversed SDEs we simply replace k with — n in the above 
equations, and generate new, uncorrected noises. 

Figure[2illustrates the time-reversal test for the calculation 
of the dynamics of the quantum average of the X -quadrature. 
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The initial coherent state has n — (a) a) = 100, and each 
stochastic trajectory was evolved forward in time until tr = 
0.5, using a constant diffusion gauge of g = 1.6. The Hamil- 
tonian was then negated and the system evolved again for the 
same period. We observe a reversal in the X-quadrature dy- 
namics back to its initial value to within statistical error, even 
though each stochastic trajectory has evolved with uncorre- 
cted noises at every point in the time-evolution. Of course, 
there are random features in each trajectory that are not time- 
reversed, but these change the distribution in a way that does 
not affect observables. This remarkable property is due to the 
overcompleteness of the quantum mechanical basis of coher- 
ent states, which permits the same physical state to be repre- 
sented in more than one way in terms of coherent states. 

While certainly not small, a one-hundred dimensional 
Hilbert space is accessible with current computers. We have 
therefore repeated this calculation for a much larger mean bo- 
son number equal to Avogadro's number n — 6.022 x 10 23 
— a truly macroscopic number of particles. This requires us 
to make use of more sophisticated gauge methods, with 



g(r, a) = 




(9) 

Further details on choices of gauges will be published else- 
where, also see 1 24 ] . 

Here the total Hilbert space dimension is astronomically 
large, and well beyond the capacity of any known or planned 
digital computer. The dynamical evolution obtained from the 
analytic result is a Gaussian amplitude decay, quite different 
to the usual exponential decay of a damped system. While the 
amplitude decay appears to correspond to information loss, in 
fact there is information stored in quantum correlations, which 
can be recovered through time-reversal. 

In this situation the physical collapse time is very short in- 
deed — orders of magnitude less than the revival time — and 
there are many oscillations of the X -quadrature. We therefore 
perform the calculation in a rotating frame, and the envelope 
of the oscillations is plotted in Fig.^b). The time reversal is 
implemented at tr = 0.5, and again we can see the revival 
of the initial state. While this situation is somewhat idealised, 
it demonstrates a fully quantum calculation for a macroscopic 
particle number. 

The important result to note is that in both cases the time- 
reversed quadrature mean agrees with the initial value to 
within the sampling error-bars. The error-bars can be reduced 
by including more trajectories in the calculation. 

For these calculations the time-reversal test illustrates a 
powerful yet counterintuitive feature of stochastic simulations 
— they can be useful for simulating unitary (reversible) quan- 
tum dynamics, despite the irreversible nature of stochastic 
processes. The examples demonstrate a stochastic simulation 
of a quantum revival, a uniquely unitary feature. As discussed 
previously, the quadrature variables will display a true revival 
to their initial values at r rcv = \ffi. By time-reversing the cal- 
culation after the initial collapse we induce the revival early. 




FIG. 2: (Color online) A representation of the broadening of the 
gauge distribution function for n — 100. The plots are histograms 
of the quantity log 10 \aQ\ for the stochastic trajectories at a number 
of points in time, (a) The time-reversed case with the Hamiltonian 
negated at tr = 0.5 and further evolved to r = 2tr. (b) The for- 
ward time evolution to r = 2tr. The initial distribution at r = for 
both graphs is a delta function, however it broadens with time to span 
several orders of magnitude at r = 2tr, despite being another repre- 
sentation of the initial quantum state in (a). Note that the logarithm 
of \aQ\ is binned. 

However, the stochastic trajectories themselves continue to 
diffuse under time-reversal — they do not simply retrace their 
forward-time path. Hence, although it seems natural to asso- 
ciate irreversible stochastic process with irreversible quantum 
dynamics (such as those of an open quantum system), clearly 
this intuition is unnecessarily restrictive. 

The calculation can equivalently be discussed in terms of 
the distribution function. During both the forwards and back- 
wards time dynamics the gauge distribution G(a) evolves ac- 
cording to a Fokker-Planck equation with positive-definite dif- 
fusion and therefore can only broaden in phase space. This is 
clearly illustrated in Fig. [2] The gauge distribution function 
that is recovered after the reversal in Fig. [2 a) is not the same 
as the initial condition, but is still equivalent to the original 
density operator. This final distribution is less compact than 
the original, but nonetheless will have identical moments cor- 
responding to the normally-ordered operator averages of the 
initial state. For this example we have sampled a gauge dis- 
tribution G(a) at 2tr that is equivalent to the initially chosen 
delta function G(a) = S 2 (a- P*)S 2 (fl~l)S 2 (a- y/n)- This 
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is precisely what must happen at the true revival time r = 
Fig. |2jb) shows the behaviour of the distribution function for 
the forward time evolution for comparison. 

Although the calculations are presented as a method of test- 
ing quantum simulations, an early revival could be observed 
experimentally with a BEC in an optical lattice using the 
phenomenon of a Feshbach resonance lEH This allows 
the tuning of both the magnitude and sign of the interaction 
strength in atomic Bose gases, represented by our parameter 
k, using an applied magnetic field. As the setup of Greiner 
et al. 12311 uses only optical potentials for the observation of 
revivals, an early revival experiment could be performed sim- 
ply with the addition of a precisely-controlled homogeneous 
magnetic field. 

We note that while the calculation presented is for a quan- 
tum phase space method, the time-reversal test is applicable to 
any quantum simulation technique. It is not a sufficient test in 
itself, since a time-reversible simulation could have other sys- 
tematic errors. However, it has the great advantage that time- 
reversibility is an exact property of unitary quantum dynamics 
even when no other exact properties are known. We believe 
that the current status of calculating the dynamics of quantum 
many-body systems is similar to the situation in the early days 
of studying classically chaotic systems on a computer. By def- 
inition such systems display sensitive dependence on initial 
conditions, and so it is difficult to estimate the errors in the 
calculated dynamics 1 19]. It was partly due to time-reversal 
tests that such calculations became convincing. 

In summary, we have presented a simple yet powerful test 
for unitary quantum dynamics, and demonstrated its use to 
verify the results of a stochastic numerical simulation in a 
macroscopically large Hilbert space. Such tests are crucial 
for these demanding calculations. An increasing variety of 
quantum dynamical techniques are now becoming available, 
and it is important to have reliable tests of their accuracy — 
especially since no analytic solutions exist in many cases of 
interest. 

We gratefully acknowledge research support from the Aus- 
tralian Research Council. 
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